【网络情报高阶实战案例】如何使用开源卫星数据进行情报项目调查

自从DigitalGlobe和Planet Labs等公司提供的高分辨率数据出现以来，卫星报道的数量如雨后春笋般涌现。尽管分辨率较低，但开源数据依旧是一种有效且及时的数据来源，还没有得到充分利用。

今天，这篇文章分享一些有关访问、理解和处理开源卫星数据的基础知识。并引导完成一些基本的案例示例——让初学者逐步进阶更高级的开源卫星数据技术的应用。        

1、了解卫星图像的工作原理

不同的卫星向地球发送不同的图像。区别包括分辨率（图像的清晰度）、它们产生的波段的数量和类型以及它们的更新频率。

解决事项

在分辨率、波段能力和可用性之间找到平衡点

图片来源：Mark Corcoran，牛津大学路透社研究所奖学金论文       

多久可以用一次？

时间因卫星而异，主要是单个卫星进入轨道和重访特定区域所需的时间。

图片来源：Mark Corcoran，

路透社研究所奖学金论文，牛津大学         

什么是光谱带？

卫星所能采样的电磁波谱中辐射的波段数(可见光、红外线、紫外线、微波、X射线等)

谱带决定了可以对数据进行何种类型的分析；

图片来源：Mark Corcoran，路透社研究所奖学金论文，牛津大学

         以哨兵二号卫星图像数据为例，从可见光、近红外到短波红外，共有十三个不同的光谱波段，其中10米处有四个光谱带，20米处有六个光谱带，60米处有三个光谱带。

将这些波段视为双筒望远镜的一种形式，可以发现原本会隐藏在数据中的事物。这些频段的正确组合是关键。可以在数据上运行各种脚本（可以使用不同的波段组合）（在本地计算机或Sentinel hub上）。        

2、教程：需要准备哪些？

1. Python 3.6

2. 一个合适的Tif阅读器（如果想下载光栅文件，则需要Tif阅读器）

3. Jupyter 笔记本和各种python包

4. Sentinel hub的免费帐户（在python教程中找到使用介绍）

使用 Sentinel Hub 浏览器工具搜索

对于部分技术专家，可能会因为使用浏览器应用程序的想法而有点反感。

对于探索和调查，EO 浏览器是一个不错的选择（如果想更进一步的话，“ Sentinel Playground ”的卫星数量较少，但提供了一种更容易探索的方式）。         

开源卫星平台可能会提供在工作流中使用python的有限选项。Sentinel Hub 在这方面运行了一些有用的选项。         

以下是EO浏览器提供的数据列表以及使用它们的理由：

卫星名称	应用描述
哨兵一号	海陆监测应急响应气候变化
哨兵2	土地覆盖图、土地变化检测图、植被监测、烧毁地区监测
哨兵3	地表地形观测海洋和陆地表面颜色观测和监测。Sentinel-3 OLCI仪器确保 Envisat Meris的连续性
哨兵5P	监测空气中一氧化碳（CO）、二氧化氮（NO2）和臭氧（O3）的浓度。监测紫外线气溶胶指数 (AER_AI) 和云的各种地球物理参数 (CLOUD)
Landsat5/7和8 的ESA档案	植被监测土地利用土地覆盖图变化监测 Landsat 8 的全球覆盖范围 - Envisat Meris和旧数据
Proba-V	土地覆盖植被生长的观察气候影响评估水资源管理农业监测和粮食安全估计内陆水资源监测和跟踪沙漠和森林砍伐的稳定蔓延。
MODIS	在全球范围内监测陆地云层海洋颜色（ESA）
GIBS	NASA 提供了包含600多颗卫星的全球图像浏览器服务

3、初级案例——追踪野火火灾情况

2018年，创纪录的火焰在美国加利福尼亚州肆虐时，发现并报告了野火的突然扩散和破坏。这可能不是最后一次。专家称，此类火灾在不久的将来可能会再次出现。

免费提供的首选资源构成了Landsat 8数据（来源美国地质调查局提供）以及Sentinel-2数据。

Sentinel-2 在可见光和红外光谱部分提供比其开源同事更高分辨率的图像 ，非常适合监测植被、土壤和水覆盖、内陆水道和沿海地区的任务。        

实战演练步骤

1. 转到EO浏览器—注册并登录（免费）

2. 选择Sentinel-2

3. 通过将云覆盖范围限制在 30% 来缩小数据收集范围。

4. 发现美国加利福尼亚州的野火，该野火在 2018年7月至 8月期间达到顶峰（它们在全州范围内的影响非常全面，您应该不难发现云层）         

2018年火灾的可能例子：

Natchez火灾（2018年7月20日）：41.956°N 123.551°W Carr火灾（2018 年7月28日）：40.6543°N 122.6236°W Mendocino Complex 火灾（2018年7月29日）：39.243283°N 123.103367°W弗格森火灾（7月14, 2018): 37.652°N 119.881°W                   

接下来，要渲染一个特定的波段组合，以更清楚地看到地面上位置上发生的情况。        

复制“野火脚本”：        

// Visualizing (wild)fires in Sentinel-2 imagery
// For use in Sinergise EO Browser (http://apps.sentinel-hub.com/eo-browser)
// https://pierre-markuse.net/2018/04/30/visualizing-wildfires-burn-scars-sentinel-hub-eo-browser/
// Pierre Markuse (@pierre_markuse)
// Wildfire and burn scar visualization in Sentinel-2 images V2.0.0
// Twitter: Pierre Markuse (@pierre_markuse)
// CC BY 4.0 International - https://creativecommons.org/licenses/by/4.0/

function a(a, b) {return a + b};
function stretch(val, min, max) {return (val - min) / (max - min);}
function satEnh(rgbArr) {
    var avg = rgbArr.reduce((a, b) => a + b, 0) / rgbArr.length;
    return rgbArr.map(a => avg * (1 - saturation) + a * saturation); }
function highlightBurnscar(val, oLow, oHigh, deSat, darken) {
    if ((B12 + B11 > 0.05) && (val > 0)) {
        if (((B8A - B12) / (B8A + B12)) > oLow) {
            saturation = saturation - deSat;
            stretchMax = stretchMax + darken;
        } else {
            if (((B8A - B12) / (B8A + B12)) <= oHigh) {
                noFire[0] = noFire[0] + 0.2 * val;
                noFire[1] = noFire[1] + 0.05 * val;
            } else {
                noFire[0] = noFire[0] + 0.15 * val;
                noFire[1] = noFire[1] + 0.15 * val;
            }
        }
    }
}
function indexMap(ind, lVal, mVal, hVal, cont, dir, pal) {
  var col1=GREEN;var col2=YELLOW;var col3=RED;  
  if (pal == 1) {col1=CBL;col2=CBM;col3=CBH;} 
  if (pal == 2) {col1=OWNL;col2=OWNM;col3=OWNH;}         
    var lValCol = col1;
    var mValCol = col2;var hValCol = col3;
    if (dir == 1){
      lValCol = col3;hValCol = col1;
    }
    if (cont == 0) {
        if (ind <= lVal) return lValCol; if ((ind > lVal) && (ind < hVal)) return mValCol; if (ind >= hVal) return hValCol;
      } else {
        return colorBlend(ind, [lVal, mVal,hVal], [lValCol,mValCol,hValCol]);
    }
}
function blend(bArr1, bArr2, opa1, opa2) {
    return bArr1.map(function(num, index) {
        return (num / 100 * opa1 + bArr2[index] / 100 * opa2);
    });
}
function applyEnh(bArr) {
    highlightBurnscar(burnscarHighlight, burnscarThresholdLow, burnscarThresholdHigh, burnscarDesaturateBackdrop, burnscarDarkenBackdrop);
    return satEnh([stretch(bArr[0], stretchMin, stretchMax), stretch(bArr[1], stretchMin, stretchMax), stretch(bArr[2], stretchMin, stretchMax)]);
}
var BLACK = [0.0, 0.0, 0.0];
var RED = [0.9, 0.1, 0.1];
var YELLOW = [0.9, 0.9, 0.1];
var GREEN = [0.0, 0.6, 0.0];
var CBL  = [0/255, 80/255, 0/255];
var CBM  = [120/255, 120/255, 230/255];
var CBH  = [70/255, 195/255, 255/255];
var OWNL = [0.0, 0.0, 0.0];
var OWNM = [0.0, 0.0, 0.0];
var OWNH = [0.0, 0.0, 0.0];
// Visualization style of the different fire zones
var Fire1OVL = [stretch((2.1 * B04 + 0.5 * B12), 0.01, 0.99) + 1.1, stretch((2.2 * B03 + 0.5 * B08), 0.01, 0.99), stretch(2.1 * B02, 0.01, 0.99)];
var Fire2OVL = [stretch((2.1 * B04 + 0.5 * B12), 0.01, 0.99) + 1.1, stretch((2.2 * B03 + 0.5 * B08), 0.01, 0.99) + 0.25, stretch(2.1 * B02, 0.01, 0.99)];
var Fire3OVL = [stretch((2.1 * B04 + 0.5 * B12), 0.01, 0.99) + 1.1, stretch((2.2 * B03 + 0.5 * B08), 0.01, 0.99) + 0.5, stretch(2.1 * B02, 0.01, 0.99)];
// Band combinations (To get quicker processing you should comment out all those you are not using in the Settings further down)
var NaturalColors = [2.9 * B04, 3.1 * B03, 3.0 * B02];
// var EnhancedNaturalColors = [2.8 * B04 + 0.1 * B05, 2.8 * B03 + 0.15 * B08, 2.8 * B02];
// var NaturalNIRSWIRMix = [2.1 * B04 + 0.5 * B12, 2.2 * B03 + 0.5 * B08, 3.0 * B02];
// var NIRSWIRColors1 = [2.6 * B12, 1.9 * B08, 2.7 * B02];
var NIRSWIRColors2 = [2.4 * B12, 1.7 * B8A, 2.2 * B05];
// var NIRSWIRColors3 = [0.5 * (B12 + B11) / 4 / B07, 0.8 * B8A, 1 * B07];
// var NIRSWIRColors4 = [2.0 * B12, 1.1 * B11, 1.6 * B08];
// var FalseColor = [B08 * 2, B04 * 2, B03 * 2];
// var NatFalseColor = [B12 * 2.6, B11 * 2, B04 * 2.7];
// var Vegetation = [B11 * 2.4, B8A * 2, B04 * 2.9];
// var PanBand = [B08, B08, B08];
// var NBR8A12 = indexMap((B8A - B12) / (B8A + B12), -0.8, -0.4, 0.0, 1, 1, 1);
// var NDVI = indexMap((B08 - B04) / (B08 + B04), -0.4, -0.2, 0.0, 1, 1, 1);
// Settings
// Fire (hot spot) visualization
var fire1 = Fire1OVL;
var fire2 = Fire2OVL;
var fire3 = Fire3OVL;
// Used band combinations and mixing
var layer1 = NIRSWIRColors2;
var layer2 = NaturalColors;
var layer1Amount = 0;
var layer2Amount = 100;
// Influence contrast and saturation
var stretchMin = 0.00;
var stretchMax = 1.00;
var saturation = 1.00;
// Fire sensitivity (Default = 1.00), higher values increase fire (hot spot) detection and false positives
var fireSensitivity = 1.00;
// Burn scar visualization
var burnscarHighlight = 0.00;
var burnscarThresholdLow = -0.25;
var burnscarThresholdHigh = -0.38;
var burnscarDesaturateBackdrop = 0.25;
var burnscarDarkenBackdrop = 0.25;
// Manually influence RGB output
var manualCorrection = [0.00, 0.00, 0.00];
// Image generation and output
noFire = blend(layer1, layer2, layer1Amount, layer2Amount);
finalRGB = applyEnh(noFire).map(function(num, index) {
    return num + manualCorrection[index];});
return (a(B12, B11) > (1.0 / fireSensitivity)) ?
    (a(B12, B11) > (2.0 / fireSensitivity)) ? fire3 :
    (a(B12, B11) > (1.5 / fireSensitivity)) ? fire2 : fire1 :
   finalRGB;

定位加州野火（2018年8月）         

另外一个有趣的角度：可看出消防员如何成功地控制/隔离火势         

Pierre Markuse的例子：

消防飞机引起的阻燃剂https://flic.kr/p/Wt8Vzo         

Pierre Markuse的代码脚本：

// VERSION=3
// QuickFire V1.0.0 by Pierre Markuse (https://twitter.com/Pierre_Markuse)
// Made for use in the Sentinel Hub EO Browser (https://apps.sentinel-hub.com/eo-browser/?)
// CC BY 4.0 International (https://creativecommons.org/licenses/by/4.0/)

function setup() {
  return {
    input: ["B01","B02","B03","B04","B08","B8A","B11","B12","CLP", "dataMask"],
    output: { bands: 4 }
  };
}

function stretch(val, min, max) {return (val - min) / (max - min);} 

function satEnh(arr, s) {
   var avg = arr.reduce((a, b) => a + b, 0) / arr.length;
   return arr.map(a => avg * (1 - s) + a * s);
}

 function layerBlend(lay1, lay2, lay3, op1, op2, op3) {
    return lay1.map(function(num, index) {
     return (num / 100 * op1 + (lay2[index] / 100 * op2) + (lay3[index] / 100 * op3));
    });
  }  

function evaluatePixel(sample) {
  const hsThreshold = [2.0, 1.5, 1.25, 1.0];
  const hotspot = 1;
  const style = 1;
  const hsSensitivity = 1.0;
  const boost = 1;

  const cloudAvoidance = 1;
  const cloudAvoidanceThreshold = 245;
  const avoidanceHelper = 0.8;

  const offset = -0.000;
  const saturation = 1.10;
  const brightness = 1.00;
  const sMin = 0.01;
  const sMax = 0.99;

  const showBurnscars = 0;
  const burnscarThreshold = -0.25;
  const burnscarStrength = 0.3;

  const NDWI = (sample.B03-sample.B08)/(sample.B03+sample.B08);
  const NDVI = (sample.B08-sample.B04)/(sample.B08+sample.B04);
  const waterHighlight = 0;
  const waterBoost = 2.0;
  const NDVI_threshold = -0.15;
  const NDWI_threshold = 0.15;
  const waterHelper = 0.2;

  const Black = [0, 0, 0];
  const NBRindex = (sample.B08-sample.B12) / (sample.B08+sample.B12); 
  const naturalColorsCC = [Math.sqrt(brightness * sample.B04 + offset), Math.sqrt(brightness * sample.B03 + offset), Math.sqrt(brightness * sample.B02 + offset)];
  const naturalColors = [(2.5 * brightness * sample.B04 + offset), (2.5 * brightness * sample.B03 + offset), (2.5 * brightness * sample.B02 + offset)];
  const URBAN = [Math.sqrt(brightness * sample.B12 * 1.2 + offset), Math.sqrt(brightness * sample.B11 * 1.4 + offset), Math.sqrt(brightness * sample.B04 + offset)];
  const SWIR = [Math.sqrt(brightness * sample.B12 + offset), Math.sqrt(brightness * sample.B8A + offset), Math.sqrt(brightness * sample.B04 + offset)];
  const NIRblue = colorBlend(sample.B08, [0, 0.25, 1], [[0/255, 0/255, 0/255],[0/255, 100/255, 175/255],[150/255, 230/255, 255/255]]);
  const classicFalse = [sample.B08 * brightness, sample.B04 * brightness, sample.B03 * brightness];
  const NIR = [sample.B08 * brightness, sample.B08 * brightness, sample.B08 * brightness];
  const atmoPen = [sample.B12 * brightness, sample.B11 * brightness, sample.B08 * brightness];
  var enhNaturalColors = [0, 0, 0];
  for (let i = 0; i < 3; i += 1) { enhNaturalColors[i] = (brightness * ((naturalColors[i] + naturalColorsCC[i]) / 2) + (URBAN[i] / 10)); }

  const manualCorrection = [0.00, 0.00, 0.00];

  var Viz = layerBlend(URBAN, naturalColors, naturalColorsCC, 10, 40, 50); // Choose visualization(s) and opacity here

  if (waterHighlight) {
    if ((NDVI < NDVI_threshold) && (NDWI > NDWI_threshold) && (sample.B04 < waterHelper)) {
     Viz[1] = Viz[1] * 1.2 * waterBoost + 0.1;
     Viz[2] = Viz[2] * 1.5 * waterBoost + 0.2;
    }
  } 

  Viz = satEnh(Viz, saturation);
  for (let i = 0; i < 3; i += 1) {
    Viz[i] = stretch(Viz[i], sMin, sMax); 
    Viz[i] += manualCorrection[i];  
  }

  if (hotspot) {  
    if ((!cloudAvoidance) || ((sample.CLP<cloudAvoidanceThreshold) && (sample.B02<avoidanceHelper))) {
     switch (style) {
       case 1:
        if ((sample.B12 + sample.B11) > (hsThreshold[0] / hsSensitivity)) return [((boost * 0.50 * sample.B12)+Viz[0]), ((boost * 0.50 * sample.B11)+Viz[1]), Viz[2], sample.dataMask]; 
        if ((sample.B12 + sample.B11) > (hsThreshold[1] / hsSensitivity)) return [((boost * 0.50 * sample.B12)+Viz[0]), ((boost * 0.20 * sample.B11)+Viz[1]), Viz[2], sample.dataMask]; 
        if ((sample.B12 + sample.B11) > (hsThreshold[2] / hsSensitivity)) return [((boost * 0.50 * sample.B12)+Viz[0]), ((boost * 0.10 * sample.B11)+Viz[1]), Viz[2], sample.dataMask];  
        if ((sample.B12 + sample.B11) > (hsThreshold[3] / hsSensitivity)) return [((boost * 0.50 * sample.B12)+Viz[0]), ((boost * 0.00 * sample.B11)+Viz[1]), Viz[2], sample.dataMask]; 
       break;
       case 2:
        if ((sample.B12 + sample.B11) > (hsThreshold[3] / hsSensitivity)) return [1, 0, 0, sample.dataMask]; 
       break;
       case 3:
        if ((sample.B12 + sample.B11) > (hsThreshold[3] / hsSensitivity)) return [1, 1, 0, sample.dataMask]; 
       break;
       case 4:  
        if ((sample.B12 + sample.B11) > (hsThreshold[3] / hsSensitivity)) return [Viz[0] + 0.2, Viz[1] - 0.2, Viz[2] - 0.2, sample.dataMask];
       break;
       default:
      }
    }
  }

  if (showBurnscars) {
   if (NBRindex<burnscarThreshold) {
     Viz[0] = Viz[0] + burnscarStrength;
     Viz[1] = Viz[1] + burnscarStrength;
   }
  }

  return [Viz[0], Viz[1], Viz[2], sample.dataMask];
}                 

如果在指定的时间范围内成功发现了一场野火，应该会发现黄红色斑点。

重要的是：不要将这些解释为火焰。尽管显示了它，但可能看到的不是真正的火灾，而只是红外覆盖——在某种程度上，这与活跃的火灾和热点一致。      

进阶练习

当时显示新的野火蔓延到肯尼亚山（经度，纬度：-0.152739，37.309095）周围

你可以尝试应用“野火脚本”去调查肯尼亚山野火情况。        

MSI

其他波段组合可用于说明存在野火风险的潜在区域。植被干燥度就是此类指标之一。水分压力指数 - 或 MSI - 可以揭示这些干燥区域 - 并有助于进行所谓的“火灾危险条件分析”。

该指数相对于其他水植被指数倒置。该值越高，水分胁迫程度越高（水分含量越少）。

// Simple Ratio 1600/820 Moisture Stress Index (abbrv. MSI)
// General formula: 1600nm / 820nm
// URL https://www.indexdatabase.de/db/si-single.php?sensor_id=96&rsindex_id=48

let index = B11 / B08;
return[index]

试一试，按照相同的步骤使用不同的波段脚本，看看你能找回什么。        

MSI脚本

//
// Simple Ratio 1600/820 Moisture Stress Index (abbrv. MSI)
//
// General formula: 1600nm / 820nm
//
// URL https://www.indexdatabase.de/db/si-single.php?sensor_id=96&rsindex_id=48
//

let index = B11 / B08;
let min = 0.058;
let max = 17.145;

// colorBlend will return a color when the index is between min and max and white when it is less than min.
// To see black when it is more than max, uncomment the last line of colorBlend.
// The min/max values were computed automatically and may be poorly specified, feel free to change them to tweak the displayed range.

let underflow_color = [1, 1, 1];
let low_color = [208/255, 88/255, 126/255];
let high_color = [241/255, 234/255, 200/255];
let overflow_color = [0, 0, 0];

return colorBlend(index, [min, min, max],
[
  underflow_color,
  low_color,
  high_color,
  //overflow_color // uncomment to see overflows
]);         

使用Python：

为了使用Sentinel Hub服务，注册一个Sentinel Hub帐户（此免费注册地址https://www.sentinel-hub.com/ ）。

登录到Sentinel Hub配置器。具有实例ID（长度为36位字母数字代码）的配置将已经存在。对于本教程，建议创建一个新配置（通过“添加新配置”）并将配置设置为基于“ Python 脚本模板”。

记下配置的实例 ID 并将其粘贴到INSTANCE_ID变量声明中： 

INSTANCE_ID = 'your ID from sentinel hub'  # In case you put instance ID into configuration file you can leave this unchanged

%reload_ext autoreload
%autoreload 2
%matplotlib inline
import datetime
import numpy as np

import matplotlib.pyplot as plt
from sentinelhub import WmsRequest, WcsRequest, MimeType, CRS, BBox
def plot_image(image, factor=1):
    """
    Utility function for plotting RGB images.
    """
    fig = plt.subplots(nrows=1, ncols=1, figsize=(15, 7))

    if np.issubdtype(image.dtype, np.floating):
        plt.imshow(np.minimum(image * factor, 1))
    else:
        plt.imshow(image)

from sentinelhub import CustomUrlParam


# Tip: if you want to insert the coordinates from google, you will need to set
# the first two coordinates for the upper left corner (-122.41, 39.31)
# and second two (-122.75, 39.55) will refer to lower right corner of the box
# Lastly: lat long from Google maps needs to be switched around (e.g. for lower corner 
# google maps will give you '39.55, -122.75'; you need to switch out around to read -122.75, 39.55)

betsiboka_coords_wgs84 = [-122.41, 39.31, -122.75, 39.55]
betsiboka_bbox2 = BBox(bbox=betsiboka_coords_wgs84, crs=CRS.WGS84)


my_url = 'https://raw.githubusercontent.com/sentinel-hub/custom-scripts/master/sentinel-2/markuse_fire/script.js'

evalscripturl_wms_request = WmsRequest(layer='TRUE-COLOR-S2-L1C', # Layer parameter can be any existing layer
                                       bbox=betsiboka_bbox2,
                                       time='2018-08-28',
                                       width=512,
                                       instance_id=INSTANCE_ID,
                                       custom_url_params={CustomUrlParam.EVALSCRIPTURL: my_url})

evalscripturl_wms_data = evalscripturl_wms_request.get_data()
plot_image(evalscripturl_wms_data[0])

所有请求都需要一个边界框作为 sentinelhub.geometry.BBox 的实例，并带有相应的坐标参考系统 (sentinelhub.geometry.CRS)。我们将使用 WGS84，我们可以使用 sentinelhub.geometry.CRS 中预定义的 WGS84 坐标参考系统。

现在只需提供JS evalscript的URL地址（此专用页面上提供了许多其他现成的脚本，地址：https://github.com/sentinel-hub/custom-scripts）。

让我们再次选择fire脚本（野火脚本）并提供其URL作为参数

CustomUrlParam.EVALSCRIPTURL的值。

Python 输出：

使用提供的 wildfire JS evalscript 下载 Sentinel-2 图像         

NBR

另一个专门用于检测火灾引起的生态系统严重程度的是NBR——标准化燃烧率的缩写（NBR详细介绍文档：https://un-spider.org/advisory-support/recommended-practices/recommended-practice-burn-severity/in-detail/normalized-burn-ratio）。

火灾强度代表有机物质在燃烧过程中释放的能量 (Keeley, 2009)。还指火势活跃时的火势强度。另一方面，烧伤严重程度描述了火灾强度如何影响被烧毁地区生态系统的功能。观察到的影响通常在区域内和不同生态系统之间有所不同（Keeley，2009 年）。烧伤严重程度也可以描述为一个区域被火灾改变或破坏的程度。图1显示了火灾强度和燃烧严重程度之间的差异。

进阶练习——使用NBR脚本找到被烧毁的植被

脚本代码：

// Normalized Difference NIR/SWIR Normalized Burn Ratio (abbrv. NBR)
// General formula: (NIR - SWIR) / (NIR + SWIR)
// URL https://www.indexdatabase.de/db/si-single.php?sensor_id=96&rsindex_id=53

let index = (B08 - B12) / (B08 + B12);
return[index]

Python 输出：

使用提供的 NBR JS evalscript 下载 Sentinel-2 图像             

中级案例——变得更具调查性        

依靠野火脚本找到准将Suheil al-Hassan的藏身处——“叙利亚最臭名昭著的军阀之一”（大西洋邮报）

这位有着可怕绰号“老虎”的将军Qawat Al-Nimr aka Tiger Forces的掌舵手，是叙利亚阿拉伯军队的一支精锐编队，在叙利亚内战中主要充当进攻部队。

据侵犯文件中心(VDC) 称，在叙利亚收复古塔东部的行动中，老虎部队最近执行的行动造成至少600名平民死亡，其中至少100名是儿童。

为了找出 Suheil al-Hassan 在2016年的藏身之处，执行了一项典型的情报工作，并从观看以下视频开始：

通过视频中可以看到向左飘来的烟味。在另一个序列中，看到了将军的藏身处。

Suheil al-Hassan 准将（又名老虎）的藏身处（左）

附近阿勒颇热电厂冒出的烟云         

从视频中，得知它是为“巴尔米拉与伊斯兰国的战斗”拍摄的。这可能是阿勒颇热电厂。

发电厂的谷歌地球图片显示了大火肆虐后的破坏程度（右侧烧毁的圆圈）        

一个简单的网络搜索为我们提供了一定程度的说明：

• 在 Google 上搜索“阿勒颇热电厂”。维基百科链接为我们提供了The Thermal Power Plant的 long/lat。

• 接下来，转到Google地球或 Google 地图并输入找到的坐标：' 36°10'30″ N 37°26'22″ E '。将看到在工厂右侧的一组被烧毁的塔楼。

在EO Brower上，打开道路并输入Google地图结果中的经度和纬度 (36.175000, 37.439444)（进入EO浏览器的搜索窗口）。在案例中，我们对2016年2月16日 (2016-02-16) 比较感兴趣，为此我们目睹了奇妙的烟雾。

烟雾向左移动         

接下来，像以前一样继续并应用markuse_fire脚本在Sentinel-2图像中可视化火灾（挑战：如果有信心，可以在你的Python环境中执行此操作，或者在EO浏览器窗口中模拟）。                 

现在我们有了更多的了解，我们可以推断出火灾发生时老虎的藏身之处。

谷歌地图上的哈桑藏身处（左）和阿勒颇热电厂（右）         

为了证实我们的怀疑，我们可以查看谷歌地图卫星图像，并了解到藏身处已经被炸毁。

Benjamin Strick是一位冲突、安全、武器和数字取证领域的开源专家，也是Bellingcat的调查员（他也提出了这个例子）解释说，它确实有助于显示当时哪些塔着火了。Al-Hassan在工厂中的后期图像可以证实：那天那四座塔着火了。

从太空中发现地方的具体细节有其优点。一个是在人权领域。最近的一项调查表明，卫星图像有助于从太空揭示一些地方的人权奴隶制。英国诺丁汉大学权利实验室数据项目主任 Doreen Boyd 估计，从太空中可以看到三分之一的人权奴隶制——无论是窑炉或非法矿山的伤痕，还是短暂的鱼类加工营地（可以说，高分辨率商业图像可能更适合这种调查）。

高阶案例——运行算法进行水库水位监测

关于目前一些非洲国家水资源水资源稀缺导致的紧张和战斗越来越有可能发生，我们此案例就进行水库水位资源情况监测。

使用Sentinel-2多光谱和多时相图像，编写了一个Jupyter notebook* 来检测水体的水位。

我们将在python中运行水检测算法，并在给定的时间间隔内提取单个水库的地表水位。        

这是我们需要做的：

• 定义一些水体的几何形状

• 准备和执行水检测的完整工作流程：使用SentinelHub服务下载 Sentinel-2 数据（真彩色和 NDWI 索引）并使用s2cloudless云检测器进行云检测（链接地址：https://github.com/sentinel-hub/sentinel2-cloud-detector），最后检测水。

• 可视化一段时间内的水体和水位。

• 过滤掉多云场景以改善结果。

需要准备的文件：       

• `eo-learn` — https://github.com/sentinel-hub/eo-learn

• `Water Observatory Backend` —https://github.com/sentinel-hub/water-observatory-backend

基本终端/文件设置：

$ means terminal code

-------------------

1. Create a new working directory and enter it

$ mkdir workshop
$ cd workshop

2. Check python version
python version should optimally be >3.6, if not, you have to install it (using sudo-apt or brew etc..)

$ python --version

3. Create python virtual environment
this creates an empty python virtual environment

$ python -m venv water_detection

load it by running

$ source water_detection/bin/activate

you should now have a specific instance of python running

4. Downloading and installing requirements
download eo-learn and water-observatory-backend from github

$ git clone https://github.com/sentinel-hub/eo-learn.git
$ git clone https://github.com/sentinel-hub/water-observatory-backend.git

5. install eo-learn

$ cd eo-learn
$ python install_all.py '-e'

go back

$ cd ..

6. upgrade pip, install ipython,jupyter notebook and other things

$ pip install --upgrade pip
$ pip install jupyter
$ pip install recordclass

downgrade tornado (current issues with latest release)

$ pip install tornado==5.1.1

7. download workshop notebook

$ wget https://raw.githubusercontent.com/mlubej/water_detection_notebook/master/water_level_extraction.ipynb

8. run the jupyter notebook and execute it

$ path_to_workdir/water_detection/bin/jupyter notebook

说明：与前面的示例一样：需要一个Sentinel Hub帐户。可以在Sentinel Hub网页上创建免费试用帐户。设置帐户后，登录Sentinel Hub Configurator。默认情况下，已经拥有带有实例ID（长度为 36 的字母数字代码）的默认配置。对于本教程，建议创建一个新配置 ( "Add new configuration") 并将配置设置为基于Python脚本模板。这样的配置将已经包含这些示例中使用的所有层。准备好配置后，请将配置的按照配置说明将实例 ID添加到sentinelhub包的配置文件中（配置说明链接：https://sentinelhub-py.readthedocs.io/en/latest/configure.html）。         

通过加载以下Python库来设置Python工作环境。确保按照上面的说明运行Python虚拟环境：

# set the autoreload and the inline plotting for matplotlib
%reload_ext autoreload
%autoreload 2
%matplotlib inline

# data manipulation
import numpy as np
import matplotlib
import matplotlib.pyplot as plt
from mpl_toolkits.axes_grid1 import make_axes_locatable

# image manipulations
from skimage.filters import threshold_otsu, sobel
from skimage.morphology import erosion, dilation, opening, closing, white_tophat, disk

# GIS related
import geopandas as gpd
from shapely.wkt import loads
from shapely.geometry import shape, MultiPolygon, Polygon

# eo-learn related
from eolearn.core import EOTask, EOPatch, LinearWorkflow, Dependency, FeatureType, LoadFromDisk, SaveToDisk
from eolearn.io import S2L1CWCSInput 
from eolearn.mask import AddCloudMaskTask, AddValidDataMaskTask, get_s2_pixel_cloud_detector
from eolearn.features import SimpleFilterTask
from eolearn.geometry import VectorToRaster

# Sentinel Hub
from sentinelhub import BBox, CRS

# water observatory backend
import sys
sys.path.append('./water-observatory-backend/src')
#from visualisation import plot_water_body
from geom_utils import get_bbox
from s2_water_extraction import get_water_level_optical
from visualisation import draw_multi, draw_poly

# other
import urllib.request as request
import json
from datetime import datetime
from shapely.wkt import loads

获取水体的几何图形

让我们以南非的 Theewaterskloof 大坝为例，这是一个重要的水库，为开普敦 400 万居民中的大部分人提供宝贵的资源。它是西开普省供水系统中最大的水坝，在干旱期间可能会出现低水位。有迹象表明人们对水资源短缺的意识有所增强。如何涵盖这样的主题显示了这个例子。

南非 Theewaterskloof 水库

对于Theewaterskloof Dam——或地球上任何其他大型水体——可以通过BlueDot Water Observatory API 轻松获取几何图形。

通过搜索特定水体，您可以复制IDURL 中的数字以访问相应水体的标称几何形状（即38538url 中的数字https://water.blue-dot-observatory.com/38538/2019-02-05）

BLUEDOT界面        

下载几何图形的Python代码：

ID = 38538

# function for obtaining the nominal water geometry from the water observatory API
def get_nominal_geometry(ID):
    wb_url = f'https://water.blue-dot-observatory.com/api/waterbodies/{ID}/index.html'
    with request.urlopen(wb_url) as url:
        wb_data = json.loads(url.read().decode())
        nominal_outline = shape(wb_data['nominal_outline']['geometry'])
    return nominal_outline

# utility function for plotting the geometry
def plot_geometry(geom, ax = None, **kwargs):
    
    if geom is None:
        return
    if geom.exterior is None:
        return
    
    x,y = geom.exterior.xy

    if ax is None:
        fig = plt.figure(figsize=(20,10))
        ax = fig.add_subplot(111)
        
    ax.plot(x, y, **kwargs)
    
# get the nominal geometry
geom = get_nominal_geometry(ID)

# and plot it
plot_geometry(geom)

现在我们需要这个几何体的边界框，以便下载Sentinel-2数据。定义了一个边界框并稍微膨胀它以构建一个与Sentinel Hub服务一起使用的BBox对象。BBox类也接受坐标系 (CRS)，我们在几何中使用相同的坐标系（即WGS84）。

# create BBox instance
bbox = get_bbox(geom, inflate_bbox=0.1)

# plot the BBox and the geometry
fig = plt.figure(figsize=(20,10))
ax = fig.add_subplot(111)
plot_geometry(bbox.geometry, ax)
plot_geometry(geom, ax)    

绘制 BBox 和几何图形         

准备/执行水检测的完整工作流程

Sentinel Hub 服务eo-learn. 它是一个用于 Python 机器学习的开源地球观测处理框架，它提供无缝访问和处理任何卫星舰队获取的时空图像序列的能力。

eo-learn作为一个工作流工作——一个工作流由一个或多个任务组成。每个任务在称为 EOPatch 的区域的一小块区域上完成特定的工作（下载数据、计算波段组合等）。EOPatch 是 EO 和非 EO 数据的容器。

让我们定义一个工作流来下载和获取水检测所需的数据。我们将下载 RGB 波段，以便实际可视化水体的真彩色图像。此外，我们将下载将用于水检测的NDWI波段组合（归一化差异水指数(https://custom-scripts.sentinel-hub.com/sentinel-2/ndwi/)）。它被定义为

归一化差水指数公式        

其中B3和B8分别是绿色和近红外Sentinel-2波段。                 

Next: 一些将在工作流中使用的自定义任务的定义代码脚本         

# calculate fraction of pixels with non-zero values
def coverage(array):
    return 1.0 - np.count_nonzero(array)/np.size(array)

# a function to return valid data for the area as a union of pixels with non-zero values and pixels that contain no clouds
class ValidDataPredicate:
    def __call__(self, eopatch):     
        return np.logical_and(eopatch.mask['IS_DATA'].astype(np.bool), 
                              np.logical_not(eopatch.mask['CLM'].astype(np.bool)))

# definition of a task to calculate and add the valid coverage scalar to the EOPatch
class AddValidDataCoverage(EOTask):
    def execute(self, eopatch):
        
        vld = eopatch.get_feature(FeatureType.MASK, 'VALID_DATA')
        
        cvrg = np.apply_along_axis(coverage, 1, np.reshape(vld, (vld.shape[0], vld.shape[1]*vld.shape[2])))
        
        eopatch.add_feature(FeatureType.SCALAR, 'COVERAGE', cvrg[:,np.newaxis])
        return eopatch

# definition of a task for water mask and water level detection
class WaterDetector(EOTask):
    def execute(self, eopatch):
        
        results = [get_water_level_optical(date, eopatch.data['NDWI'][idx,...,0], geom, bbox, simplify=True) for idx, date in enumerate(eopatch.timestamp)]
        
        df = list([x['geometry'] for x in results])
        gdf = gpd.GeoDataFrame(geometry = df, crs = {'init': eopatch.bbox.crs.ogc_string()})
        gdf['TIMESTAMP'] = eopatch.timestamp
        
        eopatch.add_feature(FeatureType.VECTOR, 'WATER_OUTLINE', gdf)
        eopatch.add_feature(FeatureType.SCALAR, 'WATER_LEVEL', np.array([x['water_level'] for x in results])[..., np.newaxis])
        
        return eopatch

EOTasks的初始化：

# TASK for downloading RGB bands
# `TRUE-COLOR-S2-L1C` is the name of the layer defined in the Sentinel Hub configurator.
# the arguments are the resolution of the image, max cloud coverage of the whole Satellite tile, and the instance ID for your Sentinel Hub account
input_task = S2L1CWCSInput(layer='TRUE-COLOR-S2-L1C', 
                           resx='20m', resy='20m', maxcc=0.5, 
                           instance_id=None)

# TASK for downloading the NDWI band combination
# other parameters are copied from the previous task
add_ndwi = S2L1CWCSInput('NDWI')

# TASK for cloud detection
# cloud probability map (CLP) and cloud mask (CLM) are calculated at 160 m resolution in order to speed up the process
cloud_classifier = get_s2_pixel_cloud_detector(average_over=2, dilation_size=1, all_bands=False)
cloud_det = AddCloudMaskTask(cloud_classifier, 'BANDS-S2CLOUDLESS', cm_size_y='160m', cm_size_x='160m', 
                             cmask_feature='CLM', cprobs_feature='CLP', instance_id=None)

# TASK for adding a raster mask of the nominal water extent (NOMINAL_WATER) 
# raster shape is provided by an existing feature inside of the EOPatch
gdf = gpd.GeoDataFrame(crs={'init':'epsg:4326'}, geometry=[geom])
add_nominal_water = VectorToRaster(feature=(FeatureType.MASK_TIMELESS, 'NOMINAL_WATER'), 
                                   vector_data=gdf, raster_value=1, 
                                   raster_shape=(FeatureType.MASK, 'IS_DATA'), 
                                   raster_dtype=np.uint8)

# TASK for adding valid data mask to the EOPatch (mask type)
add_valmask = AddValidDataMaskTask(predicate=ValidDataPredicate())

# TASK for adding valid data coverage to the EOPatch (scalar type)
add_coverage = AddValidDataCoverage()

# TASK for water detection
water_det = WaterDetector()

输出：完成加载模型，总共使用了 170 次迭代

# initialize the workflow
workflow = LinearWorkflow(input_task, 
                          add_ndwi, 
                          cloud_det, 
                          add_nominal_water, 
                          add_valmask, 
                          add_coverage,
                          water_det)

%%time

# time interval definition
time_interval = ['2016-01-01','2019-3-1']

# execute the workflow
result = workflow.execute({
    input_task: {
        'bbox': bbox,
        'time_interval': time_interval
    },
})

# result is in the form of a dictionary
eopatch = list(result.values())[-1]                          
                          

输出：CPU 时间：用户 3 分钟 9 秒，系统：14.7 秒，总计：3 分钟 24 秒     
时间：3 分钟 23 秒         

`EOPatch` 的结构

通过键入检查结构

输入：eopatch

现在通过可视化给定时间序列中所选水体的前几张真彩色图像。在下面看到一些图像包含云，这会导致正确的水位检测出现问题。       

绘制NDWI以查看水探测器如何描绘水体轮廓：

# get aspect ratio of image for better plotting
image_ar = eopatch.mask_timeless['NOMINAL_WATER'].shape[0] / eopatch.mask_timeless['NOMINAL_WATER'].shape[1]

# plot the NDWI at different dates
fig = plt.figure(figsize=(20,15*image_ar))

for i in range(12):
    ax = plt.subplot(3,4,i+1)
    ax.imshow(eopatch.data['NDWI'][i].squeeze(), vmin = 0, vmax = 1)
    ax.axis('off')
    
plt.tight_layout(pad=0)       

归一化差水指数的绘制         

绘制带有检测到的水轮廓的真彩色图像：

def plot_waterbody(img, date, dam_poly, dam_bbox, water_extent, water_level, 
                    color_nominal='white', color_current='xkcd:lime', ax = None):

    shape = img.shape[0:2]
    
    if ax is None:
        fig, ax = plt.subplots(nrows=1, ncols=1, figsize=(10,10))
        
    ax.imshow(img,extent=[dam_bbox.min_x,dam_bbox.max_x,dam_bbox.min_y,dam_bbox.max_y])
        
    if isinstance(dam_poly, Polygon):
        draw_poly(ax,dam_poly, color=color_nominal)
    elif isinstance(dam_poly, MultiPolygon):
        draw_multi(ax,dam_poly, color=color_nominal)
        
    if isinstance(water_extent, Polygon):
        draw_poly(ax,water_extent, color=color_current)
    elif isinstance(water_extent, MultiPolygon):
        draw_multi(ax, water_extent, color=color_current)


# get aspect ratio of image for better plotting
image_ar = eopatch.data['TRUE-COLOR-S2-L1C'][0].shape[0] / eopatch.data['TRUE-COLOR-S2-L1C'][0].shape[1]

fig = plt.figure(figsize=(20,15*image_ar))

for i in range(12):
    ax = plt.subplot(3,4,i+1)
    plot_waterbody(eopatch.data['TRUE-COLOR-S2-L1C'][i], eopatch.timestamp[i], geom, bbox, eopatch.vector['WATER_OUTLINE']['geometry'][i], eopatch.scalar['WATER_LEVEL'][i], ax=ax)
    ax.axis('off')
    
plt.tight_layout(pad=0)        

真实水位与Theewaterskloof Dam大坝的轮廓进行比较，一清二楚         

绘制检测到的水位

def plot_water_levels(eopatch, max_coverage=1.0):
    fig, ax = plt.subplots(figsize=(20,7))

    dates = np.asarray(eopatch.timestamp)
    ax.plot(dates[eopatch.scalar['COVERAGE'][...,0]<max_coverage],
            eopatch.scalar['WATER_LEVEL'][eopatch.scalar['COVERAGE'][...,0]<max_coverage],
            'bo-',alpha=0.7, label='Water Level')
    ax.plot(dates[eopatch.scalar['COVERAGE'][...,0]<max_coverage],
            eopatch.scalar['COVERAGE'][eopatch.scalar['COVERAGE'][...,0]<max_coverage],
            '--',color='gray',alpha=0.7, label='Cloud Coverage')
    ax.set_ylim(0.0,1.1)
    ax.set_xlabel('Date')
    ax.set_ylabel('Water Level')
    ax.set_title('Detected Water Level')
    ax.grid(axis='y')
    ax.legend(loc='best')
    return ax

# plot the water level with no cloudy scene filtering (accept all clouds)
ax = plot_water_levels(eopatch, 1.0);       

会看到由于云的干扰，数据出现了很多波动（以灰色绘制了云覆盖范围。它与水位异常值共享相同的日期）。

现在让我们为最大云覆盖率设置2%的阈值，并过滤掉与多云场景对应的日期。通过过滤掉值大于0.02的日期来完成的

eopatch.scalar['COVERAGE']。

plot_water_levels(eopatch, 0.02);

最终通过数据显示，水位在2018年年中创下三年来的历史新低。

小结：

本文除开有地理定位的应用，更多的是需要借助更专业的技术工具，通过各种代码脚本来实现调查任务，本文仅供参考了解，用于网络情报调查练习使用。

如果这篇文章对你有所帮助，麻烦点“在看”或分享给更多的网络情报爱好者。

——END——